#include "hllc.hpp"
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>

using Catch::Approx;
using namespace hms;
using namespace hms::hllc;
using s = hllc::WaveRegion;

TEST_CASE( "swe::wavespeed [sqrt(g*d)]" ){
	REQUIRE( wavespeed(1/phys::g) == Approx(1) );
	REQUIRE( wavespeed( Array2s{1, 4}/phys::g ).isApprox( Array2s{1,2} ) );
}

TEST_CASE( "normal velocity in the star region" ){
	using namespace hms::hllc;
	REQUIRE( vNormalStar(0, 0, 0, 0) == 0 );
	REQUIRE( vNormalStar(1, 1, 1, 1) == Approx(1) );
	REQUIRE( vNormalStar( 2,  3,  5,  6) == Approx( 1.5) );
	REQUIRE( vNormalStar( 3,  2,  6,  5) == Approx( 3.5) );
	REQUIRE( vNormalStar( 2., -1., 2., 1.) == Approx(1.5) );
	REQUIRE( vNormalStar(-2.,  1., 2., 1.) == Approx(0.5) );
	REQUIRE( vNormalStar(
		Array2s{0.5, 0.5},
		Array2s{-0.5, 0},
		Array2s{1, 2},
		Array2s{2, 1}
	).isApprox(Array2s{-1, 1.25}) );
}

TEST_CASE( "shallow water wave speed in the star region (cwStar)" ){
	using namespace hms::hllc;
	REQUIRE( cwStar( 2., -1., 3., 1.) == Approx(2.75) );
	REQUIRE( cwStar(-2.,  1., 3., 1.) == Approx(1.25) );
	REQUIRE( cwStar(
		Array2s{ 2,-2},
		Array2s{-1, 1},
		Array2s{1, 3},
		Array2s{3, 1}
	).isApprox(Array2s{2.75,1.25}) );
}

TEST_CASE("HLL flux"){
	using namespace hms::hllc;
	/*
	       SR*fL - SL*fR + SL*SR*(qR - qL)
	  f* = -------------------------------
	                   SR - SL
	 *
	 * with hllFlux = f*(SL, SR, fL, fR, qL, qR).
	 * fL, fR and qL, qR can be vectors. Then f* is a vector of the same size.
	 */
	SECTION("only scalar types"){
		REQUIRE( hllFlux( 5, 7, 4, 3, 8, 9 ) == Approx(24) );
		REQUIRE( hllFlux( 7, 5, 3, 4, 9, 8 ) == Approx(24) ); // reverse
		REQUIRE( hllFlux(-2., 3., 2., 3., 4., 5.) == Approx(1.2) );
	}
	SECTION("2x1 vectors for f and q, scalar for S"){
		REQUIRE( hllFlux(
			-2., 3., Array2s{2,-2}, Array2s{-3,3}, Array2s{4.,-4}, Array2s{-5,5}
		).isApprox( Array2s{10.8, -10.8}) );
	}
	SECTION("2x2 matrices for f and q, 1x2 vectors for S"){
		ArrayNNs<2,2> fL, fR, qL, qR;
		fL <<
			 2,-3,
			-2, 3;
		fR <<
			-3, 4,
			 3,-4;
		qL <<
			 4,-.5,
			-4, .5;
		qR <<
			-5, 4,
			 5,-4;
		RowArray2s SL {-2, 0.5}, SR {3, 1};
		REQUIRE( hllFlux(SL, SR, fL, fR, qL, qR).isApprox(
			Array4s{10.8,-10.8, -5.5, 5.5}.reshaped(2,2)
		) );
	}
}

TEST_CASE("Flux terms"){
	using namespace hms::hllc;

	SECTION("pressure flux"){
		REQUIRE( pressureFlux(2) == Approx(19.62) );
		REQUIRE( pressureFlux(Array2s{2, 3}).isApprox(Array2s{19.62, 44.145}) );
	}
	SECTION("mass flux"){
		REQUIRE( massFlux(2, 3) == Approx( 6) );
		REQUIRE( massFlux(2,-3) == Approx(-6) );
		REQUIRE( massFlux(Array2s{2,1.5},Array2s{3,2}).isApprox(Array2s{6,3}) );
	}
	SECTION("momentum flux, normal"){
		REQUIRE( momentumFlux_normal(2, 3) == Approx(18 + pressureFlux(2)) );
		REQUIRE( momentumFlux_normal(2,-3) == Approx(18 + pressureFlux(2)) );
		REQUIRE( momentumFlux_normal(Array2s{2,1.5},Array2s{3,2})
			.isApprox( Array2s{18,6}+pressureFlux(Array2s{2,1.5}) )
		);
	}
}
